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SUMMARY 

We study the high-contrast biharmonic plate equation with HCT and Morley discretizations. We 
construct a preconditioner that is robust with respect to contrast size and mesh size simultaneously 
based on the preconditioner proposed by Aksoylu et al. (2008, Comput. Vis. Sci. 11, pp. 319-331). By 
extending the devised singular perturbation analysis from linear finite element discretization to the 
above discretizations, we prove and numerically demonstrate the robustness of the preconditioner. 
Therefore, we accomplish a desirable preconditioning design goal by using the same family of 
preconditioners to solve elliptic family of PDEs with varying discretizations. We also present a 
strategy on how to generalize the proposed preconditioner to cover high-contrast elliptic PDEs of order 
2k, k > 2. Moreover, we prove a fundamental qualitative property of solution of the high-contrast 
biharmonic plate equation. Namely, the solution over the highly-bending island becomes a linear 
polynomial asymptotically. The effectiveness of our preconditioner is largely due to the integration of 
this qualitative understanding of the underlying PDE into its construction. 

KEY WORDS: Biharmonic equation, plate equation, fourth order elliptic PDE, Schur complement, low- 
rank perturbation, singular perturbation analysis, high-contrast coefficients, discontinuous coefficients, 
heterogeneity. 



1. INTRODUCTION 

We study the construction of robust preconditioners for the high-contrast biharmonic plate 
equation (also referred as the biharmonic equation). The aim is to achieve robustness with 
respect to the contrast size and the mesh size simultaneously, which we call as m- and h- 
robustness, respectively. In the case of a high-contrast diffusion equation, we studied the 
family of preconditioners Backs by proving and numerically demonstrating that the same 
family used for finite element discretization [4] can also be used for conservative finite volume 
discretizations with minimal modification [6]. In this article, we extend the applicability of 
Backs even further and show that the very same preconditioner can be used for a wider 
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family of elliptic PDEs. The broadness of the apphcability of Backs has been achieved by 
singular perturbation analysis (SPA) as it provides valuable insight into qualitative nature of 
the underlying PDE and its discretizations. In order to study the robustness of Bagks, we 
use an SPA that is similar to the one devised on the matrix entries by Aksoylu et al. [4] . SPA 
turned out to be an effective tool in analyzing certain behaviors of the discretization matrix 
K{m) such as the asymptotic rank, decoupling, low-rank perturbations (LRP) of the resulting 
submatrices. LRPs are exploited to accomplish dramatic computational savings and this is the 
main numerical linear algebra implication. 

The devised SPA is utilized to explain the properties of the submatrices related to K{m). 
In particular, SPA of highly-bending block Khh{itT'), as modulus of bending m ^ oo, has 
important implications for the behaviour of the Schur complement 5'(m) of Kutii'm) in Kim). 
Namely, 

S{m) Kll - KLHKHH{m)KHL = S^o + 0{m-^) , (1.1) 

where Sao is a LRP of Kll. The rank of the perturbation depends on the number of 
disconnected components comprising the highly-bending region. This special limiting form 
of S{m) allows us to build a robust approximation of S{m)~^ by merely using solvers for K^l 
by the help of the Sherman-Morrison- Woodbury formula. 

Preconditioning for the biharmonic equation was extensively studied in the domain 
decomposition setting [23, 31] and multigrid, BPX, and hierarchical basis settings [8, 14, 22, 17, 
26, 27]. Other solution strategies were also developed such as fast Poisson solvers [20, 21] and 
iterative methods [11]. However, there is only limited preconditioning literature available for 
discontinuous coefhcients. Marcinkowski [19] studied domain decomposition preconditioners for 
the mortar type discretization of the biharmonic equation with large jumps in the coefficients. 

The high-contrast in material properties is ubiquitous in composite materials. Hence, the 
modeling of composite materials is an immediate application of the biharmonic plate equation 
with high-contrast coefficients. Since the usage of composite materials is steadily increasing, 
the simulation and modeling of composite has become essential. We witness that the utilization 
of composites has become an industry standard. For instance, light weight composite materials 
are now being used in modern aircrafts by Airbus and Boeing. There is imminent need for 
robust preconditioning technology in the computational material science community as the 
modeling and simulation capability of composites evolve. 

In [29], the Euler-BernouUi equation with discontinuous coefficients was studied for the 
kinematics of composite beams. In the beam setting, the physical meaning of the PDE 
coefficient corresponds to the product of Young's modulus and moment of inertia [28] [p. 103], 
[29]. In the biharmonic plate equation setting, the PDE coefficient represents the plate modulus 
of bending [28] [p. 406]. Nonhomogeneous elastic plates has been considered in [18] with varying 
modulus of elasticity. 

Our model problem is limited to the biharmonic equation which captures only the isotropic 
materials. The extension of our analysis to a more generalized 4-th order PDE is widely open. 
Such PDEs have an important role in structural mechanics as they are used in modeling 
anisotropic materials. Plane deformations of anisotropic materials were studied in [24], but 
extension to simultaneously heterogeneous and anisotropic case needs to be further explored. 
Grossi [13] has studied the existence of the weak solutions of anisotropic plates. The coercivity 
of the bilinear forms has also been established which may lay the foundations for our future 
work related to LRPs. 



2009; 00:1-19 



ROBUST PRECONDITIONERS FOR HIGH-CONTRAST BIHARMONIC EQUATION 



3 



The remainder of the article is structured as follows. In §2, we present the underlying high- 
contrast biharmonic plate equation and the associated bilinear forms. Subsequently, the effects 
of high-contrast on the spectrum of stiffness matrix and its subblocks are also discussed. Since 
the proposed preconditioner is based on LRP, in §3, we study the LRP of the limiting Schur 
complement as in (1.1). In §4, we present the aforementioned SPA and reveal the asymptotic 
qualitative nature of the solution. In particular, the solution over the highly-bending region 
converges to a linear polynomial as m — s- cx). In §5, we introduce the proposed preconditioner 
and prove its effectiveness by establishing a spectral bound for the preconditioned system. In 
§6, a strategy is presented on how to generalize the proposed preconditioner to cover high- 
contrast elliptic PDEs of order 2k, fc > 2. In §7, the m- and /i-robustness of the preconditioner 
are demonstrated by numerical experiments. 



2. THE UNDERLYING PDE AND THE LINEAR SYSTEM 




Figure 1. Q — Qh U Ql where Qh and Ql are highly- and lowly-bending regions, respectively. 
We study the following high-contrast biharmonic equation for the clamped plate problem: 

(2.1) 



U = dnU 



f in nc M^, 

on on. 



We restrict the plate bending process to a binary regime (see Figure 1) in which the coefficient 
a is a piecewise constant function with the following values: 



a(x) 



m > 1, X e r^if , 

1, xer^L- 



It is quite common to idealize the discontinuous PDE coefficient a by a piecewise constant 
function [7, 16]. In the case of high-contrast diffusion equation, Aksoylu and Beyer [1] showed 
that the idealization of diffusivity by piecewise constant coefficients is meaningful by showing 
a continuous dependence of the solutions on the diffusivity; also see [2]. A similar justification 
can be extended to the high-contrast biharmonic plate equation. 

2.1. Bilinear forms for the biharmonic equation 

In the theory of elasticity, potential energy is defined by using rotationally invariant functions. 
For plates, the potential energy is given by [9, p. 30]: 



J^i'^) '■— 7^ [ [{trace iJess}^ + 2((T — 1) det ffess] dx — j fv dx, 
2 Jo Jo 



(2.2) 
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where Hess is the Hessian, 



diiv di2V 

d2lV d22V 



Hess — 

The bihnear form corresponding to energy minimization in (2.2) is given by: 

a{u,v) :— / a [V^uV^v + (I — a){2di2udi2V ~ diiud22V — d22udiiv}'\ dx, (2.3) 

where < cr < 1/2 is the Poisson's ratio. Note that the straightforward bihnear form associated 
to (2.1) is obtained by using Green's formula: 

/ (a V^u) u = / aV'^uV^v dx + / adnV^uv dj - / aV^udnV dj. (2.4) 
Jn Jn Jan Jan 

We see that both (2.3) and (2.4) contain the so-called canonical bilinear form, a(u,w), 
associated to the biharmonic equation (2.1): 



d{u,v) := / a\7^u\/^v dx. 
Jn 



(2.5) 



When u,v € Hg^fl), both bilinear forms a{u,v) and d(u,v) correspond to the strong 
formulation (2.1) due to second Green's formula and the zero contribution of the below term: 



/ (1 — (T){29i2u9i2f — 9iiu922W — d22udiiv} dx. 
Jn 

2.2. Effects of high- contrast on the spectrum 



(2.6) 
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Figure 2. The HCT discretization of the biharmonic equation with m = 10^". (Left) The spectrum of 
the stiffness matrix K. (Right) Spectrum of the diagonally scaled stiffness matrix. Notice the 3 small 
eigenvalues of order 0{m~^) corresponding to the kernel of the Neumann matrix, span{]^^, x^, y^}. 
The plot of the two of smallest eigenvalues overlap because they are roughly of the same magnitude. 
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Roughness of PDE coefficients causes loss of robustness of preconditioners. This is mainly 
due to clusters of eigenvalues with varying magnitude. Although diagonal scaling has no effect 
on the asymptotic behaviour of the condition number, it leads to an improved clustering in 
the spectrum. The spectrum of diagonally scaled stiffness matrix, A, is bounded from above 
and below except three eigenvalues in the case of a single isolated highly-bending island. On 
the other hand, the spectrum of K contains eigenvalues approaching infinity with cardinality 
depending on the number of DOF contained within highly-bending island. For the case of 
HCT discretization with m = lO^'^, we depict the spectra of K and A and their subblocks in 
Figure 2. Clustering provided by diagonal scaling can be advantageous for faster convergence 
of Krylov subspace solvers especially when deflation methods designed for small eigenvalues 
are used; for further discussion see [5]. 

Utilizing the matrix entry based analysis by Graham and Hagger [12] for linear FE, in [6], 
the authors extended the spectral analysis to cell-centered FV discretization and obtained 
an identical spectral result for A. Namely, the number of small eigenvalues of A depends 
on the number of isolated islands comprising the highly-bending region. We observe a similar 
behaviour for the biharmonic plate equation where the only difference is that for each island we 
observe three small eigenvalues rather than one. The three dimensional kernel of the Neumann 
matrix is responsible for that difference; see §3. A similar matrix entry based analysis can be 
applied to discretizations of the plate equation, but this analysis is more involved for HCT and 
Morley discretizations than that for linear FE. Hence, we exclude it from scope of this article. 



3. DISCRETIZATIONS AND LOW-RANK PERTURBATIONS 

We consider an _ff^-conformal and also an i?^-nonconformal Galerkin finite element 
discretization; Hsieh-Clough- Tocher (HCT) [10] and Morley [25] elements, respectively. Let 
the linear system arising from the discretization be denoted by: 

K{m)x^b. (3.1) 

is decomposed with respect to magnitude of the coefficient value as 

Sl^H/fUr^L, (3.2) 

where Q.h and 0,^ denote the highly- and lowly-bending regions, respectively. DOF that lie on 
the interface, F :— f2/f nfi^, between the two regions are included in ^h- When m-dependence 
is explicitly stated and the discretization system (3.1) is decomposed with respect to (3.2), 
i.e., the magnitude of the coefficient values, we arrive at the following 2x2 block system: 

KHH{m) Khl 
_ Klh Kll 

There are important properties associated to the Khh block in (3.3): It is the only block that 
has m-dependence, and furthermore, a matrix with low-rank kernel can be extracted from it. 
Our preconditioner construction is based on LRPs from this extraction. Next, we explain how 
to extract the so-called Neumann matrix and why a(u, v) is the suitable bilinear form for that 
purpose. 

By rewriting (2.3) as the following 

a{u,v) — / a [aV^uV^v + {I — (7){diiudiiv + d22ud22V + 2di2udi2v}] dx, (3.4) 





XH 




bH 




XL 




bL 



(3.3) 
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we see that 

a{v,v) = acr |jV2i;|li^(t^) +a(l - cr)|w|^2(f^) 

> (3.5) 

The inequahty (3.5) has important imphcations. Namely, a{v,v) is Vp^ (f2)-coercive where 
Vpiifl) C H'^{n) is a closed subspace such that V'p^{il) n Pi = and Vi denotes the set 
of polynomials of degree at most 1. Furthermore, (3.5) immediately implies that a{v,v) is 
iJQ(ri)-coercive. 

Let T'' be the triangulation of ft. Based on T'*, we define the associated discrete space 
Vj]; (fi) such that Vji n Vi = 0. A precise definition of the Khjj block in the stiffness matrix 
in (3.1) is given by: 

:-a(^J„V^^), 

where (t)^,^'^ € V^{^h) C Hq{^Ih) are the basis functions. We define the Neumann matrix 
Nhh as follows: 

where (f>']:j,ip'^ S Vp_^{il.H)- Since a(-, •) is (ri)-coercive, this implies by (3.5) that 

keiAfHH = 'PiIuh = span{l^,x^,y^}. (3.6) 
Hence, KHHi'm) has the following decomposition: 

KHH{m)^mAfHH + R, (3.7) 

where R is the coupling matrix corresponding to DOF on the interface F. Now, we are in a 
position to reveal the resulting main numerical linear algebra implication. As to — » oo, the 
limiting Schur complement 6*00 in (1.1) becomes a rank-3 perturbation of K^l. This result 
relies on the fact that the inverse of the limiting Khh is of rank-3; see (4.1). This is due to 
the fact that Mhh has a rank 3 kernel whose (normalized) discretization is given by: 

eif := [lij,^ff,y„]. (3.8) 



4. MAIN SINGULAR PERTURBATION ANALYSIS RESULTS 

Lemma 4.1. The asymptotic behaviour of the submatrices in (5.1) is given by the following: 

KHH{m)-^ = eH?7~'e^ + 0(m-i), (4.1) 

S{m) = KLL~{KLLeH)'n-\e'HKLL) + 0{m-^), (4.2) 

KLHKHH{m)-^ = {KLLeH)r]-^e'H + 0{m-^), (4.3) 

where 

T]:=e*HKHHeH. (4.4) 

Proof Since Mhh is symmetric positive semidefinite, using (3.6) we have the following spectral 
decomposition where nn denotes the cardinality of DOF in CIh'- 

Z'NhhZ = diag(Ai, . . . , A„^„3, 0, 0, 0), (4.5) 
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Z'KhhMZ 



where en is defined in (3.8). Using 



where {A^ : i = 1, . . . ^uh} is a non-increasing sequence of eigenvalues of J\fHH and Z is 
orthogonal. Since, the eigenvectors corresponding to the zero eigenvalues are discretization of 

the polynomials l,x, and y, we can write Z — Z \ 

(3.7), we have: 

'm dia.g{Xi, . . . ,XnH-3) + Z*RZ Z^Ren 
e^fjRZ e^fRen 

A(m) 6 



5' 7]} 

To find the limiting form of KHni'm)^^ note that 
A(m) = m diag(Ai, . . . , A„^_3) + 

= m diag(Ai, . . . , A„^_3) (^1 + m^^ diag(Aj;\ . . . , A,7^_3)Z*i?Z^ 

Then, 

\\A{mr% < max.<„,_3 A-^^ 

1-m 1 maxi<„^_3 A^ \\Z'^RZ\\2 
for sufficiently large m, we can conclude the following: 

We proceed with the following inversion: 



A(to) S 
(5* rj 



Uirn) V{m) U{mY 



where 

U{m) 
V{m) 

Then, (4.7) implies that 



0* 1 
A(to)"^ 







0* 



77 - ^*A(to) 



U{m) 
V{m) 

Combining the above results, we arrive at 



O 

0* r/-i 



" A(m) 5 ' 


-1 


' " 


5* r] 




0* rr^ 



+ 0(m~i). 



and, by (4.6), we have 



= Z 



Z' + 0{m-^) 



O 

0* T]-^ 



(4.6) 



(4.7) 



(4.8) 



2009; 00:1-19 



8 



BURAK AKSOYLU AND ZUHAL YETER 



KHH{m) 



(4.9) 



which proves (4.1) of the Lemma. 

Parts (4.2) and (4.3) follow from simple substitution and using (5.2). | 

Remark 4.1. // we further decompose DOF associated with fin into a set of interior DOF 
associated with index I and interface DOF with index T, we obtain the following block 
representation of Khh^ 

Kii{m) Kir{m) 
Kri{m) Krrim) 

The entries in the block Krrin^) oltc assembled from contributions both from finite elements 

in Q,H and i-e. KYvifn) — A^^vi^) + ^it ■ 

We further write Ch in block form; en = (ej , e^)*. Finally we note that the off-diagonal 
blocks have the decomposition: 

Klh ^ [0 Klt ] = ^HL- (4.10) 
Therefore, the results of Lemma 4-1 can be rewritten as the following: 

Sim) = Kll ~ (i^irer) (ef^if^^^er)^' (e^i^ri + O(m-i), 
KLHKHHimy^ = (i^irer) (efi^^p'er)"'e^ + 0(m-i). 
4.I. Qualitative nature of the solution 

We advocate the usage of SPA because it is a very effective tool in gaining qualitative insight 
about the asymptotic behavior of the solution of the underlying PDE. Through SPA, in 
Lemma 4.1, we were able to fully reveal the asymptotic behaviour of the submatrices of K in 
(5.1). This information leads to a characterization of the limit of the underlying discretized 
inverse operator. We now prove that the solution over the highly-bending island converges to 
a linear polynomial. In other words, x'^ G span en- This is probably the most fundamental 
qualitative feature of the solution of the high-contrast biharmonic plate equation. 

Lemma 4.2. Let en as in (3.8). Then, 

XH{m)^ CH CH + 0{m-^), (4.11) 
where c/f is a 3 x 1 vector determined by the solution in the lowly-bending region. 

Proof We prove the result by providing an explicit quantification of the limiting process based 
on Lemma 4.1: 

XL{m) = 5"^(m) {6l - is:LHi^^^(m)6/f} 

= S^^{bL-KLH{eHV-^e'H)hH} + 0{m-^) 

=: xf+0{m-^), 

XH{in) = K^\j{m) {bu - KHLXL{m)} 

= eHri-^e^H{^H-KHLX^} + 0{m-^) 

=: CH CH + 0(to-1). 
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5. CONSTRUCTION OF THE PRECONDITIONER 



The exact inverse of K can be written as: 



Ihh 




III 



Khh 








Ihh 
-Kj.hK 





III 



(5.1) 



where Ihh and III denote the identity matrices of the appropriate dimension and the Schur 
complement S is exphcitly given by: 



S{m) = Kll - KLHKHHi^)I^HL- 



(5.2) 



Let the Hmit in (4.1) be denoted by K'^l^ :— enrj ^e^. Based on the above perturbation 
analysis, our proposed preconditioner is defined as follows: 



BAGKs{m) 



Ihh -K^„Khl 
III 



KHH{my 








s- 



LHH 







-KlhK^h III 



(5.3) 



We need the following auxiliary result to be used in the proof of Theorem 5.1 which 
characterizes the spectral behaviour of the preconditioned system. 

Lemma 5.1. For sufficiently large m, we have 

^HH' = effry-i/24 + 0(m-V2), (5.4) 

where rj is the 3x3 SPD matrix independent of m defined in (4.4). 

Proof We start by writing down the spectral decomposition of KHH{m) 

Q{mYKHH{m)Q{m) = diag(^i(m), . . . , /i„„_3(m), /i„^_2(m), /i„^_i(m), /^„^(m)), 

where {niim) : i — l,...,n//} denotes a non-increasing ordering of the eigenvalues of 
KHnifn). Since KHH{m) is SPD, we have Hi{m) > for all i < nn- We use the main fact 
that eigenvalues and eigenvectors of a symmetric matrix are Lipschitz continuous functions of 
the matrix entries [15, 30]. 

By (4.5) and (4.8) in Lemma 4.1, we give the following spectral decomposition: 



(5.5) 



Note that rj in (4.6) is a 3 x 3 symmetric, and hence, diagonalizable matrix. We proceed towards 
a fully diagonalized form of the limiting i4r^^(m). For that, we use the diagonalization of r]~^: 

77-1 = zh, Mffi + zh^ IJ-hI + fJ-n] ZHy ■ 

Therefore, we have the following expression for the last term in (5.5): 

enV^^en = [zh^ zh^ z//J diag(^^|, ' MhJ) [zh^ zh^ zhJ\ (5.6) 

where 

[zH, ZH^ ZHy] ■■= [eHi en^ e^/J [zh, zh^ i^J 
[eHi,eH^,eHy] := e^. 
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Now by substituting (5.6) in (5.5), we have the following spectral decomposition which 
corresponds to the fully diagonalized version: 

-^ffffC™) = ZiOzl + ... + Z„„_3 Z*,^_3 + ZH, Mffi z\j^ + ZH^ flH^ Z*H^ + ZH^ fiHy Z^Hy + C'( 

= : Zeo diag(0, . . . , 0, Mffl, Mffi, /%i) Zl + 0{m-^). 

The expression in (5.7) also implies the convergence of the eigenvectors of Khh{iti)'- 

Q(m) = Zoo + O(m-1). (5.8) 

Note that Zoo differs from Z in (4.5) only in the last three columns due to diagonalization of 
V- 

From (5.7), we obtain a characterization of the largest three eigenvalues of Kniii'm)^^'- 

f^nH-2{m)-^ = /i^| + 0(m-i) (5.9a) 
Ai„„_i(m)-i = (1^1+0(771-') (5.9b) 
/i„„(m)-i = fi^l + 0{m-^) . (5.9c) 

Using (5.7) and (5.9), we arrive at the following: 

diag(^i(m)"^/^, . . . , ^„„_3(m)"^/^, ^„„_2("^)"^^^, Ai«„-i("i)"^/^, Ain„ (m)"^/^) 
= diag(0, . . . , 0, vr^j\ix-^j^) + 0(m-i/2). (5.10) 

By using (5.10) and (5.8), we arrive at the desired result: 

A'ifff(m)-i/2 = Q(m)diag(/ii(m)-i/2,...,/i„„(m)-i/2)Q(m)* 

= Zoo diag(0, . . . , 0, yr^l\ Mh'^', A^^^^") C + 0{mr^l^) 

r 11-/ -1/2 -1/2 -1/2n r it , -i/2\ 

= [zifi Zff, Zff J diag(/i^/ ) [z_f/, z//^ +0{m, ' ) 

= en r,-i/2g*^ + 0(^^-1/2)^ 



The following theorem shows that Backs is an effective preconditioner for m 3> 1. 
Theorem 5.1. For sufficiently large m, we have 

a{BAGKs{m) K{m)) C [1 - cm-^/^ 1 + cm"^/^] 
for some constant c independent of m, and therefore 

^{,Bagks{.tti) K{m)) - 1 + 0{m-^'''). 

Proof Let us factorize the preconditioner as Backs — L^L with 



L := 



KHHim)-'/^ 
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where the Umiting Schur complement S'(m) and KLHKfj\j is denoted by Sao and 
respectively. We can easily show that 

a{BAGKsK) = a{LKL') = a{I + E). (5.11) 

Note that 

P^hKhmPlh - PlhKhl - KLH{eHV'^e%KHHeHri-^e'H - eHri-^e^H^HL - 0. (5.12) 

We give a step of the operation leading to (5.11). Using (5.12), the (2, 2)-th block entry of the 
LKL^ reads: 

The other entries of LKL* can be computed in a similar way. 
Using (5.4), we have 

Elh = S^^/^KLH{lHH-eHV'^e*HKHH)eHrr'^''e'H + 0{m-^/^) = 0{m-^/^). 

Hence p{E), the spectral radius of E, is ©(m^^^^), which together with (5.11) completes the 
proof. I 



6. GENERALIZATION TO ELLIPTIC PDES OF ORDER 2k 

In essence, the biharmonic plate equation preconditioner is an extension of the construction 
for the diffusion equation. It is possible to generalize this construction to a family of elliptic 
PDEs of order 2k, k > 2. We present how to obtain LRPs from associated bilinear forms. We 
choose a different perspective than the one in Section 3. We start with a canonical bilinear 
form and show the modification it needs to go through in order to construct LRPs. 
Let the generalized problem be stated as follows: Find u E i/p (fi) such that 

Tku := (-1)'= V'^ {ak V'u) = / in ft. (6.1) 

The straightforward bilinear form associated to (6.1) is obtained by application of Green's 
formula k times: 

/ V'' {akV''u)v dx ^ / akV''uV''v dx + boundary terms. (6.2) 
Jn Jn 

Then, we define a bilinear form corresponding to (6.1) which can be seen as a generalization 

of the canonical bilinear form in (2.5): 

dk{u,v) := / ttfc V'^uV'^t; da;. (6.3) 
Jn 

Without modification, ak{-, •) cannot lead to LRPs because dk{v, v) is not i/p (ri)-coercive. This 
is due to the fact that dk{v, v) — for v E Vk-i H HQ{il). Hence, the stiffness matrix induced 
by (6.3) has a large kernel involving elements from 'Pk_i n V'^ which indicates that extraction 
of a Neumann matrix with a low-dimensional kernel is impossible. In order to overcome this 
complication, we utilize a modified bilinear form: 

ak{u,v) = dk{u,v) + (1 - CTfe) dk{u,v). 

The bilinear form should maintain the following essential properties: 
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1. (f2)-coercive. 
2- '^Pfc_i(n)-coercive. 

3. Corresponds to a strong formulation giving TkU in (6.1) precisely, 

where Vp^._j(f2) is a closed subspace such that V-Pf^_^{^) nVk-i = and Vk-i denotes the set 
of polynomials of degree at most fc — 1. 

The above properties (1) and (2) will be immediately satisfied if the generalization of (3.5) 
holds for the modified bilinear form: 

ak{v,v)> Cklv^Hk^Qy (6.4) 

A similar construction of the Neumann matrix can be immediately generalized as follows: 

:= afc(0^„V^). 

The low-rank perturbations arise from the following decomposition of K^^jj{m): 

where 77'^'^^ :— e^"* K^^e^ . LRP is produced by e^' G T-'k^i because the rank is equal to the 
cardinality of the basis polynomials in V^_i. 

Due to (2.6), a2(-, •) in (2.3) corresponds to the strong formulation T2 exactly. Let us denote 
the strong formulation to which afe(-, •) corresponds by Tk- We have Tk = Tk, A: = 1, 2 for the 
high-contrast diffusion and biharmonic plate equations, respectively: 

ai(w,ti) := (Vu, ai Vf ) 

a2(w,f) := 0-2 (V^u, a2 V^w) H- q;2 (1 - CT2)|w|^2(n) 

However, for general fc, afc(-, •) may not correspond to T^. In addition, one may need more 
general boundary conditions if similar zero contributions in (2.6) can be obtained for general fc. 
Further research is needed to see if such boundary conditions are physical. Currently, it is also 
unclear for which applications such general PDEs can be used. However, there are interesting 
invariance theory implications when one employs bilinear forms corresponding to rotationally 
invariant functions compatible to energy definition in (2.2). This allows a generalization of the 
energy notion and may be the subject for future research. For further information, we list the 
relevant bilinear forms that are composed of rotationally invariant functions derived by the 
utilization of invariance theory. 

a3{v,v) := CT3 (V^t;,a3 V^w) -I- as (1 - cr3)|w|^3(n) 
a4{v,v) := 0-4 (V'*u,a4 V^w) -I- a4 (1 - a-4)|i'lH4(f^) -I- a4 74|V^i;|^2(n)- 
Note that the above bilinear forms satisfy (6.4). 
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7. NUMERICAL EXPERIMENTS 

The goal of the numerical experiments is to compare the performance of the two 
preconditioners: AGKS and MG. The domain is a unit square whose coarsest level triangulation 
consists of 32 triangles. We consider the case of a single highly-bending island located at the 
region [1/4,2/4] x [1/4,2/4] consisting of 2 coarsest level triangles. For an extension to the case 
of multiple disconnected islands, one can refer to [4, Sections 3 and 4]. The implementation 
of HCT and Morley discretizations are based on Pozrikidis' software provided in [28]. The 
problems sizes of HCT and Morley discretizations are 131, 451, 1667, 6403 and 81, 289, 1089, 
4225 for levels 1,2,3 and 4, respectively. 

We denote the norm of the relative residual at iteration i by rA"^^: 

where r*-*^ denotes the residual at iteration i with a stopping criterion of rr^'^ < 10~^. In Tables 
I-VIII, preconditioned conjugate gradient iteration count and the average reduction factor 
are reported for combinations of preconditioner, smoother types, and number of smoothing 
iterations. The average reduction factor of the residual is defined as: 



We enforce an iteration bound of 60. If the method seems to converge slightly beyond this 
bound, we denote it by 60^^, whereas, stalling is denoted by oo. 

We use Galerkin variational approach to construct the coarser level algebraic systems. The 
multigrid preconditioner MG is derived from the implementation by Aksoylu, Bond, and 
Hoist [3]. We employ a V(l,l)-cycle with point symmetric Gauss-Seidel (sGS) and point Gauss- 
Seidel (GS) smoothers. A direct solver is used for the coarsest level. 

By exploiting the fact that 6*00 in (1-1) is only a LRP of Kll, we can build robust 
preconditioners for Soo in (5-3) via standard multigrid preconditioners. (1.1) implies that 

•S'oo = Kll - vq'^^v^ , 

where v := Kj^hCh- If Mj^i^ denotes a standard multigrid V-cycle for Kj^i^, we can construct 
an efficient and robust preconditioner S^^ for Soo using the Sherman-Morrison- Woodbury 
formula, i.e. 

S-^ Mll + Mllv {v-v^Mllv)-^v'^Mll- (7.1) 

Note also that we can precompute and store Mj^j^v during the setup phase. This means that 
we only need to apply the multigrid V-cycle M^^^ once per iteration. Therefore, the following 
practical version of preconditioner (5.3) is used in the implementation: 



Backs '■= 



Ihh -K^hKhl 
III 



Mhh 

5-1 



Ihh 
-KlhKWh III 



(7.2) 



We construct two different multilevel hierarchies for multigrid preconditioners Mhh in (7.2) 
and Mll in (7.1) for DOF corresponding to Q,h and Q.l, respectively. For prolongation, linear 
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Table I. AGKS + HOT + sGS + smooth number 1-5-10 



N\m 


10° 




10^ 10^ 10" 10^ 


10^ 


109 


lOlo 


smooth number = 1 



T31 24, 0.485 20, 0.447 18,0.407 17,0.371 17,0.381 16,0.337 18,0.371 16,0.362 17,0.384 

451 52,0.730 38,0.650 21,0.452 13,0.286 12,0.249 12,0.256 13,0.279 12,0.253 11,0.213 

1667 60+, 0.857 60+, 0.768 33, 0.610 20, 0.426 18, 0.401 19, 0.410 21, 0.447 19, 0.420 19, 0.417 

6403 00,0.972 60+, 0.930 60+, 0.839 45,0.692 37,0.637 36,0.636 36,0.638 36,0.635 39,0.661 



smooth number = 5 

T31 24, 0.485 20, 0.447 18,0.407 17,0.371 17,0.381 16,0.337 18,0.371 16,0.362 17,0.384 

451 40,0.664 28,0.547 15,0.330 8,0.131 6,0.054 6,0.023 4,0.014 4,0.016 4,0.012 

1667 60+, 0.786 48,0.706 24,0.490 12,0.258 8,0.091 6,0.058 5,0.035 5,0.026 5,0.024 

6403 60+, 0.947 60+, 0.862 43,0.682 21,0.427 12,0.223 8,0.091 6,0.051 6,0.052 6,0.062 

smooth number = 10 

T31 24, 0.485 20,0.447 18,0.407 17,0.371 17,0.381 16,0.337 18,0.371 16,0.362 17,0.384 

451 37,0.634 26,0.528 15,0.330 8,0.131 6,0.050 6,0.017 4,0.010 3,0.004 3,0.003 

1667 60+, 0.785 43,0.680 20,0.442 12,0.213 8,0.080 6,0.030 4,0.004 4,0.002 4,0.008 

6403 60+, 0.943 60+, 0.861 38,0.653 20,0.410 10,0.177 8,0.090 5,0.028 5,0.015 5,0.023 



Table 11. AGKS + HOT -h GS -h smooth number 1-5-10 



N\m 10° 


lOi 


102 


10^ 


10" 


10^ 


10^ 


109 


lOio 








smooth 


number 


= 1 








131 24,0.485 


20, 0.447 


18,0.407 


17, 0.371 


17,0.381 


16,0.337 


18,0.371 


16,0.362 


17, 0.384 


451 57,0.749 


49, 0.720 


27,0.538 


23, 0.494 


22,0.459 


23, 0.480 


26,0.535 


24, 0.490 


25,0.517 


1667 60+, 0.918 60+, 0.880 


60+, 0.872 60+, 0.847 60+, 0.853 60+, 0.820 60+, 0.871 60+, 0.881 60+, 0.814 


6403 oo, 1.001 


00,0.991 


00,0.958 


00,0.953 


00,0.964 


00,0.971 


00,0.980 


00,0.977 


00,0.985 








smooth 


number 


= 5 








131 24,0.485 


20,0.447 


18,0.407 


17,0.371 


17,0.381 


16,0.337 


18,0.371 


16,0.362 


17,0.384 


451 37,0.644 


26,0.538 


16,0.339 


8,0.133 


6,0.053 


6,0.024 


4,0.011 


4,0.010 


4,0.014 


1667 60+, 0.786 


48, 0.706 


23,0.494 


12,0.253 


8,0.106 


6,0.060 


5,0.022 


5,0.022 


5, 0.027 


6403 60+, 0.947 60+, 0.887 


50,0.724 


22, 0.480 


12,0.253 


9,0.141 


10,0.185 


9,0.138 


10,0.163 








smooth 


number 


= 10 








131 24,0.485 


20, 0.447 


18,0.407 


17, 0.371 


17,0.381 


16,0.337 


18,0.371 


16,0.362 


17, 0.384 


451 37, 0.637 


25,0.525 


14,0.312 


8,0.131 


6, 0.050 


6,0.016 


4, 0.002 


3,0.004 


3, 0.002 


1667 60+, 0.785 


43,0.680 


20, 0.442 


12,0.213 


8,0.080 


6,0.029 


4,0.005 


4,0.002 


4, 0.006 


6403 60+, 0.946 60+, 0.861 


45,0.696 


20,0.410 


10,0.196 


8,0.085 


6,0.052 


6,0.033 


6, 0.040 
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Table III. AGKS + Morley + sGS + smooth number 1-5-10 



N\m 


10" 


IQi 


10^ 10^ lO"* 10^ 


10^ 


109 


lOlo 


smooth number = 1 



"81 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.316 10,0.189 8,0.080 6,0.030 4,0.012 4,0.004 3,0.004 3,0.001 3,0.0004 

1089 25,0.519 17,0.367 10,0.182 8,0.074 6,0.027 4,0.015 4,0.005 4,0.003 4,0.003 

4225 52,0.733 34,0.619 18,0.358 10,0.173 7,0.084 6,0.044 6,0.035 6,0.043 6,0.036 



smooth number = 5 

81 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.316 10,0.189 8,0.080 6,0.031 4,0.012 4,0.004 3,0.004 2,0.0003 3,0.0004 

1089 25,0.514 17,0. .363 10,0.181 8,0.074 6,0.024 4,0.009 4,0.001 3,0.002 3,0.001 

4225 46,0.698 27,0.546 16,0.315 10,0.152 6,0.057 6,0.018 4,0.003 4,0.002 3,0.004 

smooth number = 10 

"81 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.316 10,0.189 8,0.080 6,0.031 4,0.012 4,0.004 3,0.004 2,0.0002 3,0.0004 

1089 25,0.514 17,0.363 10,0.181 8,0.074 6,0.024 4,0.009 4,0.001 3,0.002 3,0.001 

4225 46,0.698 27,0.546 16,0.315 10,0.151 6,0.057 6,0.018 4,0.003 4,0.002 4,0.004 



Table IV. AGKS + Morley + GS + smooth number 1-5-10 



N\m 


10" 


lOi 


10^ 103 lO* lO^ 


10^ 


109 


lOlo 


smooth number = 1 



8i 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.329 10,0.189 8,0.080 6,0.031 4,0.012 4,0.006 3,0.005 3,0.003 3,0.003 

1089 28,0.550 19,0.402 10,0.192 8,0.085 6,0.043 5,0.040 5,0.037 6,0.030 5,0.039 

4225 59, 0.760 38, 0.651 20, 0.433 12, 0.249 11, 0.204 11, 0.214 11, 0.208 11, 0.210 11, 0.206 



smooth number = 5 

"81 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.316 10,0.189 8,0.080 6,0.031 4,0.012 4,0.004 3,0.004 2,0.0002 3,0.0004 

1089 25,0.514 17,0.363 10,0.181 8,0.074 6,0.024 4,0.009 4,0.001 3,0.002 3,0.001 

4225 46,0.698 27,0.0546 16,0.315 10,0.152 6,0.057 6,0.018 4,0.003 4,0.002 3,0.004 

smooth number = 10 

81 9,0.119 7,0.066 5,0.040 4,0.015 4,0.005 4,0.002 2,0.0003 2,0.0003 3,0.0002 

289 15,0.316 10,0.189 8,0.080 6,0.031 4,0.012 4,0.004 3,0.004 2,0.0001 3,0.0004 

1089 25,0.514 17,0.363 10,0.181 8,0.074 6,0.024 4,0.009 4,0.001 3,0.002 3,0.001 

4225 46,0.698 27,0.0546 16,0.315 10,0.152 6,0.057 6,0.018 4,0.003 4,0.002 3,0.003 



2009; 00:1-19 



BURAK AKSOYLU AND ZUHAL YETER 
Table V. MG + HCT + sGS + smooth number 1-5-10 



N\m 10° IQl 10^ 10* IQB 10^ lO'' 10^ 10^ 



smooth number = 1 

131 60+, 0.885 60+, 0.898 60+, 0.932 00,0.988 oo,0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 cxD, 0.963 00,0.987 ool. 014 oo, 1.050 oo, 1.086 oo, 1.106 oo, 1.172 oo, 1.081 oo, 1.091 

1667 00,0.985 oo, 1.015 oo, 1.044 oo, 1.062 oo,1.122 oo, 1.109 oo, 1.142 oo, 1.170 oo, 1.124 

6403 00,1.025 oo, 1.040 oo, 1.057 oo, 1.125 oo, 1.145 oo, 1.130 oo, 1.171 oo, 1.112 oo, 1.187 

smooth number = 5 

131 60+, 0.885 60+, 0.898 60+, 0.932 00,0.988 oo,0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 60+, 0.761 60+, 0.829 60+, 0.920 00,1.070 oo,1.084 oo, 1.120 oo, 1.174 oo, 1.118 oo, 1.166 

1667 60+, 0.854 60+, 0.923 oo, 0.999 oo, 1.038 oo, 1.0037 oo, 1.0085 oo, 1.134 oo, 1.154 oo, 1.208 

6403 60+, 0.931 oo, 0.979 oo, 0.998 oo, 1.012 oo,1.023 oo, 1.058 oo, 1.041 oo, 1.063 oo, 1.099 

smooth number = 10 

131 60+, 0.885 60+, 0.898 60+, 0.932 00,0.988 oo,0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 48,0.660 53, 0.701 60+ , 0.825 oo, 0.955 oo, 1.032 oo, 1.115 oo, 1.179 oo, 1.200 oo, 1.196 

1667 40,0.624 49, 0.680 60+, 0.797 oo, 1.001 oo,1.088 oo, 1.035 oo, 1.064 oo, 1.052 oo, 1.095 

6403 60+, 0.890 60+, 0.929 oo, 0.972 oo, 1.049 oo, 1.017 oo, 1.052 oo, 1.051 oo, 1.134 oo, 1.170 



Table VI. MG -|- HCT -h GS smooth number 1-5-10 



N\m 10° lOi 10^ 10"* 105 10^ lO'' 10* lO'' 



smooth number = 1 

131 60+, 0.885 60+, 0.898 60+, 0.932 oo, 0.988 oo, 0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 oo, 0.995 oo, 1.016 oo, 1.029 oo, 1.041 oo, 1.107 oo, 1.114 oo, 1.128 oo, 1.184 oo, 1.205 

1667 oo, 1.033 oo, 1.034 oo, 1.042 oo, 1.077 oo, 1.115 oo, 1.155 oo, 1.068 oo, 1.068 oo, 1.079 

6403 oo, 1.055 oo, 1.052 oo, 1.057 oo, 1.070 oo, 1.146 oo, 1.141 oo, 1.146 oo, 1.115 oo, 1.160 

smooth number = 5 

131 60+, 0.885 60+, 0.898 60+, 0.932 oo, 0.988 oo, 0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 59, 0.760 60+, 0.828 60+, 0.947 oo, 1.043 oo, 1.061 oo, 1.121 oo, 1.149 oo, 1.176 oo, 1.189 

1667 60+, 0.857 60+, 0.904 oo, 1.003 oo, 1.033 oo, 1.056 oo, 1.101 oo, 1.116 oo, 1.152 oo, 1.176 

6403 oo, 0.961 oo, 1.003 oo, 1.037 oo, 1.072 oo, 1.084 oo, 1.103 oo, 1.105 oo, 1.115 oo, 1.122 

smooth number = 10 

131 60+, 0.885 60+, 0.898 60+, 0.932 oo, 0.988 oo, 0.997 oo, 1.075 oo, 1.089 oo, 1.065 oo, 1.137 

451 37,0.632 42,0.667 60+, 0.778 oo, 1.060 oo, 1.081 oo, 1.103 oo, 1.161 oo, 1.194 oo, 1.200 

1667 47,0.706 57, 0.753 60+, 0.853 oo, 1.048 oo, 1.022 oo, 1.040 oo, 1.072 oo, 1.114 oo, 1.149 

6403 60+, 0.924 60+, 0.949 oo, 1.008 oo, 1.028 oo, 1.033 oo, 1.048 oo, 1.037 oo, 1.052 oo, 1.069 
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17 



N\m 10° 10^ 10^ 10* 10^ 10^ 10''' 10* 10^ 



smooth number = 1 

81 38,0.652 45,0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 3 7,0.638 45,0.679 54,0.736 oo, 1.014 oo, 1.034 oo, 1.021 oo, 1.031 oo, 1.098 oo, 0.021 

1089 50,0.724 60,0.766 60+, 0.869 oo, 1.036 oo,1.001 oo, 1.003 oo,1.002 oo, 1.116 oo,1.055 

4225 60+, 0.877 00,1.009 oo, 1.021 oo, 1.021 oo, 1.061 oo, 1.057 oo, 1.120 oo, 1.151 oo, 1.163 

smooth number = 5 

"81 38,0.652 45,0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 13,0.242 13,0.282 17,0.355 26,0.526 30,0.583 43, 0.686 60+, 0.830 oo, 1.027 oo,1.148 

1089 13,0.270 16,0.362 20,0.407 31,0.585 46, 0.691 60+, 0.817 oo, 1.025 oo,1.005 oo,1.055 

4225 34,0.621 41,0.668 51, 0.726 60+, 0.793 60+, 0.906 oo, 1.095 oo, 1.027 oo, 1.081 oo, 1.031 

smooth number = 10 

81 38,0.652 45,0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 9,0.111 10,0.135 11,0.207 14,0.307 18,0.395 23,0.465 27,0.526 43, 0.680 60+, 0.800 

1089 12,0.219 14,0.290 1 7,0.380 21,0.460 26,0.477 32,0.571 45, 0.691 60+, 0.772 oo, 0.977 

4225 31,0.593 36,0.632 44,0.682 55, 0.742 60+, 0.814 60+, 0.900 oo,1.038 oo, 1.125 oo,1.061 



Table VIII. MG + Morley + GS + smooth number 1-5-10 



N\m 10° 10^ 10^ 10" 10^ 10® lO'^ 10® 10^ 



smooth number = 1 

"81 38,0.652 45,0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 52,0.724 60+, 0.807 oo, 0.955 oo,0.989 oo, 0.982 oo, 1.043 oo, 0.990 oo, 1.027 oo,1.020 

1089 60+, 0.860 60+, 0.894 oo,0.996 oo,0.989 oo, 1.047 oo, 1.091 oo, 1.021 oo,1.036 oo, 1.173 

4225 00,0.972 oo, 1.011 oo, 1.020 oo, 1.066 oo, 1.058 oo, 1.129 oo, 1.134 oo, 1.145 oo, 1.164 

smooth number = 5 

"81 38,0.652 45, 0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 14,0.243 16,0.284 18,0.332 31,0.547 38, 0.646 60+, 0.826 oo, 1.037 oo, 1.082 oo, 1.085 

1089 16,0.364 21,0.441 27,0.517 45, 0.699 60+, 0.774 oo, 1.014 oo, 1.042 oo, 1.020 oo,1.038 

4225 39,0.652 50, 0.718 60+, 0.765 60+, 0.844 60+, 0.942 oo, 1.073 oo, 1.092 oo, 1.107 oo, 01.123 

smooth number = 10 

"81 38,0.652 45,0.694 54,0.741 60+, 0.841 60+, 0.924 60+, 0.921 oo, 1.001 oo, 1.045 oo, 0.959 

289 10,0.104 10,0.186 12,0.246 16,0.274 21,0.382 26,0.537 38,0.626 45, 0.697 60+, 0.870 

1089 15,0.289 16,0.362 20,0.405 25,0.522 30,0.519 36, 0.607 47, 0.702 60+, 0.904 oo,1.009 

4225 36,0.635 42,0.678 51, 0.721 60+, 0.830 60+, 0.919 oo, 1.027 oo, 1.048 oo, 1.003 oo, 01.114 
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interpolation is used as in [8]. The prolongation matrices Phh and Pll are extracted from 
the prolongation matrix for whole domain U, in the fashion following (3.3): 

p _ Phh Phl 
" [ Plh Pll _ ' 

As emphasized in our preceding paper [4], AGKS can be used purely as an algebraic 
preconditioner. Therefore, the standard multigrid preconditioner constraint that the coarsest 
level mesh resolves the boundary of the island is automatically eliminated. However, for a fair 
comparison, we enforce the coarsest level mesh to have that property. 

We do not observe convergence improvement when a subdomain deflation strategy based 
on the smallest eigenvalues is used as in the diffusion equation case [6]. The eigenvectors of 
the Neumann matrix, ch in (3.8), cannot approximate the eigenvectors corresponding to the 
smallest eigenvalues of Khh which are of 0(1) (see Figure 2) since the remainder matrix R in 
(3.7) is of 0(10^). Therefore, a deflation strategy utilizing en will not necessarily guarantee 
deflation of the smallest eigenvalues of K^h in the biharmonic case. 

We first observe that the Morley discretization provides faster convergence for both 
preconditioners. Then, we have the following results regarding the effect of number of 
smoothing iterations on the convergence behaviour. The convergence of MG heavily depends 
on the number of smoothing iterations, i.e., the more the smoothing iteration, the faster 
the convergence. For the HCT discretization, AGKS requires more than a single smoothing 
iteration for convergence; see Tables I and II. However, for the Morley discretization, even 
with the same minimal number of smoothing iteration, AGKS leads to convergence; see Tables 
HI and IV. The choice of 5 smoothing iterations is sufficient for AGKS to reach /i-robustness 
and its peak performance. Hence, we can conclude that AGKS clearly enjoys /i-robustness. 
In contrast, MG is not /i-robust regardless of the m value and the smoothing number; see 
Tables V, VI, VII, and VIII. MG is totally ineffective as the problem size increases for both 
discretizations, and more obviously for HCT. 

Finally, we report the m-robustness results. The loss of m-robustness of MG can be observed 
consistently for all m values; see Tables V, VI, VII, and VIII. The AGKS preconditioner 
becomes more effective with increasing m and reaches its peak performance by maintaining 
an optimal iteration count for all m > 10^. This indicates that m > 10^ corresponds to the 
asymptotic regime. Even increasing the m value from lO'^ to 10'^ reduces the iteration count 
signiflcantly, a clear sign of close proximity to the asymptotic regime. In addition, the AGKS 
outperforms MG even for m — 1. Consequently, for both discretizations, we infer that AGKS 
is TO-robust. 
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